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FOREWORD 


^  An  approach  particularly  suitable  for  phenomenological  band  structure 

calculations  of  highly  anisotropic  structures  is  developed.  The  approach  is 
based  on  a  cylindrical  multipole  expansion  along  a  high-symmetry  axis  and  the 
associated  reduction  of  the  Schrodinger  equation  into  a  set  of  coupled,  1-d 
wave  equations.  The  convergence  properties  of  the  multipole  expansion  are 
tested  for  the  low  F-L  bands  of  PbSe.  Good  convergence  is  demonstrated. 
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ABSTRACT 


A  method  for  band  structure  calculations  along  a  high-symmetry  axis  is 
introduced  in  this  report.  It  entails  a  cylindrical  coordinate  multipole 
expansion  of  the  wave  function  and  an  exact  reduction  of  the  three- 
dimensional  Schrodinger  equation  into  a  set  of  one-dimensional  wave  equations 
for  the  multipoles.  Group  theoretic  considerations  and  energy  arguments 
imply  good  convergence,  regardless  of  the  unit-cell  extension  along  the  high- 
symmetry  direction.  Band  calculations  of  a  test  case,  PbSe  along  the  [111] 
direction,  demonstrate  the  good  convergence.  This  method  is  expected  to  be 
particularly  useful  for  structures  highly  anisotropic  along  a  high- symmetry 
axis,  e.g.,  superlattices. 
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SECTION  1 
INTRODUCTION 


The  band  model  forms  the  basis  for  quantitative  understanding  of 
electronic  properties  of  crystalline  solids.  It  grew  out  of  a  need  to 
reconcile  two  opposite  characteristics:  on  the  one  hand  solids,  e.g.,  metals, 
exhibit  independent  particles  (gas)  features  such  as  the  resistivity,  while  on 
the  other  hand,  the  expected  strong  short-range  electron-electron  repulsion 
and  electron-ion  attraction  imply  a  highly  correlated  system.  This  basic 
issue  was  resolved  by  the  Landau  theory  of  normal  Fermi  liquids  and  the 
concept  of  pseudopotentials.  Landau  theory  clarified  how  a  strongly 
correlated  fermion  system,  e.g.,  electrons,  responds  to  an  external  one-body 
probe  as  if  it  is  comprised  of  independent  constituents- -i . e . ,  the 
quasiparticles.  Provided  the  probe -generated  excitations  are  near  the 
chemical  potential  (i.e.,  the  Fermi  level)  the  quasiparticles  behave  as  the 
"bare"  fermions  except  for  an  effective  mass.  The  pseudopotential  concept,  on 
the  other  hand,  describes  the  apparent  weak  interaction  between  near-Ep 
electrons- -the  valence  electrons --and  the  ionic  cores. Underlying  this 
weak  interaction  is  the  "cancellation  theorem"'^  which  states  that,  at  the 
ionic -core  region,  the  strong  valence  and  core  electrons  Coulombic  repulsion 
approximately  cancels  the  strong  valence-electrons  and  atomic  ion  attraction. 
The  net  attraction  to  the  ionic  cores  is  consequently  weak.  When  combined, 
these  two  fundamental  developments  lead  to  the  band  structure  picture:  "Real" 
valence  electrons  (up  to  the  mass)  move  independently  of  each  other  in  an 
effective,  weak,  one  body,  lattice -periodic  potential,  i.e.,  the 
pseudopotential.  This  electronic  motion  is  described  in  terms  of  the 
appropriate  energy  bands.  The  resulting  theory  yields  quantitative  agreement 
with  a  large  body  of  data  such  as  the  frequency  dependent  dielectric  constant, 
band  gaps,  electronic  transport,  etc. 

Practical  implementation  of  this  physical  picture  calls  for  the 
construction  of  the  pseudopotential  and  efficient  numerical  algorithms  to 
solve  the  associated  Schrodinger  equation.  A  large  number  of  approaches 
exist^  ranging  from  empirical  pseudopotentials  and  a  small  set  of  basis  states 
to  sophisticated,  self-consistent,  nonlocal  pseudopotentials  in  ab-initio 
schemes.^  This  spectrum  of  methods  represents  a  tradeoff  between  "simple" 
phenomenological  calculations  and  realistic  complex  codes.  For  materials  with 
a  large  number  of  atoms  in  the  unit  cell,  ab-initio  calculations  are 
impractical . 

The  advent  of  superlattice  fabrication,  i.e.,  the  growth  of  layered 
structures  with  alternate  constituents,  poses  new  challenges  to  band  structure 
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calculations.  A  major  difficulty  is  the  colossal  size  of  the  unit  cell, 
comprised  of  all  unit  cells  in  a  superlattice  period.  As  mentioned  above, 
most  band  structure  codes  cannot  handle  very  large  unit  cells.  Several 
methods  have  been  developed  to  address  this  difficulty.  A  popular 
phenomenological  approach  is  the  envelope -function  method,  supplemented  by  a 
prescription  for  the  boundary  conditions  across  the  interface. °  Another 
approach  is  ab-initio  calculations  in  a  small  supercell. ^  Other  methods,  such 
as  the  complex  energy  band  method,®  combined  features  of  both  types  of 
"limits,"  In  this  work  we  introduce  yet  another  method,  specifically  tailored 
for  phenomenological  calculations  of  highly  anisotropic  structures  along  a 
symmetry  axis.  The  method  is  predicated  on  the  usage  of  an  empirical 
pseudopotential^  and  is  limited  to  bands  along  a  high-symmetry  direction. 

Our  starting  point  is  the  observation  that  for  an  epitaxial  superlattice 
there  is  a  preferred  direction  in  space,  i.e.,  the  growth  direction,  which  is 
invariably  also  a  high  symmetry  axis.  This  feature  suggests  using  cylindrical 
coordinates  where  the  z-axis  coincides  with  the  system's  high  symmetry  axis. 
This  leads  to  the  cylindrical  coordinates  multipole  expansion  (Section  2). 
Explicitly,  any  lattice -periodic  function  f{r),  e.g.,  the  pseudopotential  or 
the  periodic  part  of  the  wave  function,  can  be  exactly  expanded  in  terms  of 
cylindrical  multipoles  f^(g,z'): 

00 

=  X  X  J^(g  (1-1) 

g 

where  p,<|),z'  are  the  standard  cylindrical  coordinates,  J^(x)  is  the  Bessel 
function  of  integer  order  i,  and  g  varies  over  a  infinite  series  of  positive 
numbers  determined  by  the  geometry  of  the  reciprocal  vector  space.  As 
Equation  (1-1)  shows,  the  4)  and  p  (transversal)  dependence  is  carried  by 
standard  functions.  Consequently,  the  analysis  is  focused  on  solving  for  the 
z' -dependent  multipoles -- the  "guts"  of  f(r).  To  keep  the  analysis  simple, 
spin  degrees  of  freedom  are  omitted  at  this  stage.  A  subsequent  paper  covers 
the  appropriate  extension. 

The  usefulness  of  expansion  Equation  (1-1)  hinges  on  its  convergence 
properties,  i.e.,  on  how  many  multipoles  are  needed  to  describe  relevant 
bands.  The  other  important  issue  is  the  simplicity  of  the  equations  governing 
the  multipoles.  The  arguments  below  and  the  detailed  discussion  in  the 
following  sections  show  that  Equation  (1-1)  has  good  convergence  properties 
and  that  the  multipoles  satisfy  a  set  of  one-dimensional  wave  equations. 

The  good  convergence  of  expansion  Equation  (1-1)  is  a  consequence  of 
severe  limitations  on  the  range  of  the  i  and  g  labels  due  to  group- theoretical 
and  energy  considerations.  To  demonstrate  these,  we  analyze  the  case  of  a 
rocksalt  structure  along  the  [111]  direction  (superlattice  or  bulk  material). 
As  shown  in  Section  3,  the  [111]  symmetry  axis  implies  that  only  two  i  values 
need  be  considered!  This  renders  Equation  (1-1),  in  effect,  to  an  expansion 
in  g.  The  g  label,  in  turn,  controls  the  number  of  transversal  nodes  through 
the  Bessel  function  oscillations.  Specifically,  for  i>0,  J^(g  p)  rises  from 
J^(0)  -  0  to  a  maximum  at  p-(f  +  3/2)/g  and,  henceforth,  oscillates  with  a 
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q 

period  1/g.  (For  i  -  0,  the  latter  oscillations  start  at  p  -  0.) 
Consequently,  for  a  given  energy  band,  all  g  values  up  to  a  cutoff  value  g^ 
contribute.  A  crude  estimate  of  g^  is  obtained  as  follows.  In  view  of  the 
above,  "low"  energy  bands  have  the  least  niimber  of  p-oscillations  in  the  unit 
cell.  To  obtain  zero  p-oscillations  within  the  unit  cell,  it  is  necessary 
that  (ijji  +  3/2)g^  »  a-j.,  where  a.p  is  a  typical  transversal  extension  of  the 
unit  cell  and  is  a  typical  £  value.  (Since  only  two  i  values  enter,  is 

well  defined.)  For  the  test  case  of  bulk  PbSe  (Section  5) 
f jjj  =  3 ,  aj  -  a/2  =  3A  where  a  is  the  lattice  constant.  This  gives 
g^  =  1.5  A~^  which  is  approximately  the  second  g  value  (Section  3,  Table  3-1). 
Thus,  for  this  example,  only  about  four  terms  in  Equation  (1-1)  need  be 
considered.  The  above  considerations  apply  equally  to  a  superlattice  in  the 
[111]  direction  since  the  transversal  projection  of  the  lattice  constant  is 
the  same  for  a  bulk  material  and  a  superlattice.  Hence,  the  good  convergence 
is  expected  to  hold  for  the  latter. 

The  convergence  properties  of  Equation  (1-1)  are  tested  numerically  for 
the  above  case  in  the  entire  [111]  segment  of  the  Brillouine  zone.  The 
results  demonstrate  good  the  convergence,  with  six  multipoles  or  less 
(Section  3) . 

The  other  result  of  this  work  is  the  derivation  of  an  exact  set  of  wave 
equations  for  the  wave  function  multipoles  (Section  4).  Explicitly 

r  2  ^"*0  -♦  d^  1 

-gj  +  — H  E(k)  +  -^U^(gQ.z')  + 

*-  ft  dz  J 

3  r  2ra_  1 

+  I  I  n,<gO>6P'gF)  -  —2  '^m^®P’^'Ti-m^SF-2'>  -  0  (^-2) 

gp.gp  ’  '■ft  ■> 

where  ^'^(g,z')  and  "^niCgp-z')  are  the  wave  function  (unknown)  and 
pseudopotential  (known)  multipoles,  respectively,  and  the  A-coef f icients 

m^gO’gp’gp^  geometrical  coefficients  which  depend  on  the  lattice  point 

^roup  and  the  symmetry  of  z'  axis.  To  our  knowledge,  the  geometrical 
A-coeff icients  are  new  constructs  which  embody  the  transversal  momenta  vector 
addition  geometry.  The  one -dimensionality  of  Equation  (1-2)  is  advantageous 
when  dealing  with  a  unit  cell  which  extends  considerably  in  the  high-symmetry 
axis  direction.  It  can  be  solved  using  any  algorithm  for  coupled  ordinary 
differential  equations. 

The  cylindrical  multipoles  method  developed  here  has  its  advantages  and 
limitations.  As  pointed  out  above,  its  primary  advantage  is  good  convergence, 
particularly  for  highly  anisotropic  unit-cells.  Thus,  for  example,  the  method 
can  handle  a  large  number  of  atoms  in  the  unit  cell.  On  the  down  side,  our 
method  is  not  valid  throughout  the  Brillouine  zone,  hence  unsuitable  for  self- 
consistent  calculations. 

The  paper  is  organized  as  follows:  In  Section  2  we  introduce  the 
multipole  expansions  for  the  pseudopotential  and  the  wave  function.  Section  3 
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is  a  discussion  of  the  symmetry  constraints  on  the  multipole  expansion.  In 
Section  4  we  derive  the  basic  set  of  coupled  wave  equations.  The  numerical 
results  for  bulk  PbSe  bands  in  the  T-A-L  direction  are  given  in  Section  5. 
Section  6  is  a  summary  and  brief  discussion. 
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SECTION  2 

THE  CYLINDRICAL  MULTIPOLE  EXPANSION 


As  in  other  versions  of  the  empirical  pseudopotential  approach,  it  is 
necessary  to  solve  the  three-dimensional  Schrodinger  equation  with  a  local 
periodic  potential.  Unlike  other  approaches,  however,  we  introduce  a 
cylindrical  coordinate  system  from  the  outset  by  choosing  the  z-axis  along  a 
high- symmetry  direction  of  the  crystal.  This  choice  has  been  advocated  in  the 
INTRODUCTION  as  useful  in  conjunction  with  superlattices,  surfaces,  and  for 
anisotropic  materials.  In  this  Section  we  start  to  develop  the  approach.  For 
the  sake  of  specificity,  it  is  hereafter  assumed  that  the  high-symmetry  axis 
is  the  [111]  direction  and  the  crystal  structure  is  fee. 


The  transformation  from  the  (x,y,z)  coordinates  to  the  primed 
coordinates,  where  z'  is  parallel  to  [111],  is  given  by 


x' 

lf2“  iFF  [2“ 

2  3  2J3  ‘J3 

X 

y' 

■= 

0 

22 

y 

Z  ' 

1_  1- 

z 

>[3"  J3~  J3" 

»  « 

or 

(?')  -  R  (r) 


(2-la) 


(2-lb) 


The  cylindrical  coordinates  corresponding  to  (x’,  y',  z')  are  conventionally 
defined  by 


x'  =  p  cos4),  y'  “  p  sin(}>,  z'  =  z*  (2-2) 

where  ^  is  the  azimuthal  angle.  It  is  this  cylindrical  coordinates  system 
that  is  used  throughout.  By  the  same  token,  denoting  the  reciprocal  lattice 
vectors  in  the  unp^ime^  coqrdinates  by  (5  -  (G^^,  G^,  G^),  the  vector's 
components  2'  -  (G^^,  G  ,  G^)  in  the  primed  coordinates  are  given  by 
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(G')  -  R  (G)  (2-3) 

since  2  transforms  as  a  polar  vector.  The  fee  reciprocal  lattice  vectors  are 
given 


a 

(  -  n^^  n2  -t-  n^ ,  np 

-  n2  +  np ,  np  n2  -  np) 

np,  n2,  np  -=  0,  ±1, 

±2,  ±3  .  .  . 

(2-4) 

and  a  is  the  lattice  constant. 


To  introduce  the  multipole  expansion  of 
consider  first  the  empirical  pseudo  potential 


a  lattice-periodical  function, 
v(r)^ ’ 


v(?)  =  X'  -2'  v'(G')  -v'(r') 

G  G' 


(2-5) 


where  v'(2')  =  v(2(2'))  and  the  summation  in  Equation  (2-5)  is  over  a 
restricted  set  of  2  vectors  of  the  form  Equation  (2-4).  The  coefficients  v(2) 
are  obtained  from  fitting  selected  data.  Now  the  exponential  in 
Equation  (2-5)  can  be  transcribed  to  the  chosen  cylindrical  coordinates 
Equation  (2-2) : 


G'*r'  -  p  gp  cos  (4i  -  &(G.p))  +  G^z'  (2-6) 

9  .  # 

where  the  transyersal  reciprocal  vector  G.j.  is  Gj  =  (G^.Gy)  and  the  azimuthal 
shift  angle  6(2^)  and  length  gp  are  defined  by: 

cos  6(G;)  -  —  .  sin  6(G;)  - 

T  gp  T  gp 

gp  -  J(G^)2  -h  (G^)2  (2-7a) 

For  the  case  gp  -  0,  only  the  G^  component  survives  and,  hence,  we  adopt  the 
convention: 


6(G,J.)  -  0  for  gp  -  0 


(2-7b) 


The  multipole  expansion  is  derived  by  employing  the  identity 


12 


iz  cos  (|) 


I 


.i  T  /  \ 

1  (z)  e 


(2-8) 
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where  J^Cz)  denotes  the  Bessel  function  of  order  i.  Applying  Equation  (2-8) 
to  the  RHS  of  Equation  (2-5)  gives 


where 


v'(r')  -  X  X  v^(gp,z')  J^(gp  p) 
gp 


.ii<t) 


(2-9a) 


,  V  ’-G,  z'  -iiS(G') 

v^(gp,z')  =  1  X  v'(G')  e  z  e  "  T" 

G' 

I Gf I -gp 


(2-9b) 


The  sununation  in  Equation  (2-9b)  is  over  all  2'  such  that  the  length  of  2^  has 
the  prescribed  value  gp.  Also,  since  J^(0)  -  0’ 


v^(gp  -  0,z')  -  &2Q  Vq(0,z') 
Equation  (2-9c)  is  consistent  with  Equation  (2-7b) 


(2-9c) 


A  similar  analysis  can  be  applied  to  a  band  wave  function  • 


.13 


T  /■♦x  ik»r  ^ 


£„£(?)  -  I  »„S  (I) 


iL*r 


nk 


(2-10) 


When  transformed  to  the  primed  coordinates,  the  lattice  periodic  function 
fnj^(r)  is  expanded  as  in  Equation  (2-9),  while  the  exponential  factor  takes 
the  form  exp[ilc*r]  -  exp[ilc'*r']  -  exp[ik'z']  for  Ic  vectors  pointing  in  the 
direction  of  the  high-symmetry  axis.  For  such  Ic  vectors- -the  only  Ic  vectors 
considered  hereafter- -  the  wave  function  multipole  expansion  is 


where 


^V(r') 

nk 


gp 


'I«^^(gp,k',z')  -  i^  X  e^^^z+  k')z'^-ii6(Lp) 


L' 

\K\  -  gp 


(2-lla) 


(2-llb) 


and  the  shift  angle  is  defined  as  in  Equation  (2-7) 
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f  0 

cos  6(L  )  =  ,  sin  &(L  )  - ^ 

T  gp  ^  Sp 

gp  -  J(L^)2  +  (1^)2  (2-12) 

As  before,  in  Equation  (2-llb)  and  when  gp  -  0  then  6(11^)  - 

i  “  0.  The  notation  in  Equation  (2-lla)  is  simplified  hereafter  by  omitting 
the  k',n  labels,  i.e.,  *ni ^ ^ ^  “  ^^(gp-^')- 

The  attractive  feature  of  a  multipole  (2-9)  and  (2-11)  is  that  it  lumps 
the  contribution  of  many  plane  waves,  all  of  which  have  the  same  transversal 
reciprocal  vector  projection  length.  Hence,  few  multipoles  may  have  the 
capacity  of  describing  well  the  wave  function.  In  addition,  the  unknown 
^^(gF,z')  depends  only  on  one  variable  and,  hence,  are  considerably  easier  to 
solve . 
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SECTION  3 

SYMMETRY  PROPERTIES 


The  multipole  expansions  Equations  (2-9)  and  (2-11)  are  general  and 
exact.  The  particular  structure  and  symmetry  of  the  crystal  determines  the 
values  of  g,  i  and  the  locus  of  the  transversal  reciprocal  vectors  2.J..  We 
discuss  these  matters  below  for  the  [111]  direction  and  fee  crystals. 

Consider  first  the  allowed  g  values  and  the  locus  of  2.^.  These  can  be 
calculated  from  Equations  (2-3)  and  (2-4),  with  result^  summarized  in 
Figure  3-1  and  Table  3-1.  Figure  3-1  shows  that  all  2.p  with  a  given  length  g 
fall  into  hexagonal  stars.  These  stars  come  either  as  singles  (Figure  3-la) 
or  as  twins  (Figure  3-lb).  In  the  latter  case  the  two  hexagons  are  , 

symmetrically  tilted  with  respect  to  by  an  angle  ®Q(g) .  Therefore,  for  all 
the  shift  angle  6(2.j.)  has  the  form 


6(G^) 


i  +  f  j)>  J  -  -2.  -1. 


+3 


(3-1) 


The  values  of  g  and  'J>Q(g)  for  the  first  few  2^  stars  are  given  in  Table  3-1. 
Note  that  the  lowest  g  values  are  sparse.  This  is  a  desirable  feature.  Since 
the  transversal  undulations  of  the  Bessel  function  J^(pg)  have  a  period 
-g'^  A,^  it  is  expected  that  the  lowest  bands  involve  g  lower  than  a 
threshold.  Therefore,  a  sparse  sequence  of  g's  implies  that  only  few  g  values 
contribute  to  the  lowest  bands. 


An  important  consequence  of  Equation  (3-1)  is  obtained  by  combining  it 
with  definition  Equation  (2-llb).  This  gives 


’^i±6^SF’^'^  -  -  e"^^®0<6F>  ♦^(gp.z')  . 

when  sin[64>Q(gp)  ]  -  0  (3-2a) 

sin[12$Q(gp) ] 

'*'i+12^®F’^’^  “  ■  '*'i+6^^F’^'^  sin[6*Q(gp)]  ’ 

when  sin[6$Q(gp)]  0  (3-2b) 

Equations  (3-2)  (and  others  not  published  here)  imply  that  the  number  of 
f - independent  multipoles  is  <6.  This  is  the  first  severe  constraint  on 
Equation  (2-11)  due  to  symmetry.  Note  also  that  the  four  lowest  g  values  in 
Table  3-1  pertain  to  Equation  (3-2a). 


3-1 


NAVSWC  TR  91-32A 


-»  o  ^ 

<  2  m 

o  <  - 

OC  Ui  5 
<-  I  O 
•-  P 
1  Q  < 

^  z  2 

V  <  o 


u.  ^ 

o  ci  r: 

UJ  Ui  Ot 

_)  Ci 
<  t  ,P 
U  t- 


U  <  z 

z  z  Q 

1/1  o  ^ 

2  o  g 

0  8  2 

Ui  tu 
i/>  Ui  JTJ 

<  I  iA 

O  K 


3-2 


FIGURE  3-1.  THE  LOCUS  OF  THE  fee  RECIPROCAL  VECTORS  PROJECTIONS 
UNTO  THE  PLANE  NORMAL  TO  THE  [111]  DIRECTION 
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TABLE  3-1.  STRUCTURE  OF  fee  g-HEXAGONS  ASSOCIATED  WITH  THE  [111]  DIRECTION 


g  -  f/a 


•C 

g(PbSe)<^) 

(A‘^)  (degrees) 

plane  wave 
index  N(n)^^^ 

Hexagon 

Multipliei 

1. 

0. 

0. 

0. 

3j 

1 

2. 

10.2606 

--IT 

1.6754 

0. 

^3j 

1 

3. 

17.7714 

- 

2.9018 

30 

3j 

2 

4. 

20.5207 

--F 

3.3507 

0. 

^3j 

1 

5. 

27.1465 

4.4326 

19.1066 

''3j 

2 

6. 

30.7813 

5.0261 

0. 

3j 

1 

7. 

35.5430 

5.8036 

30 

3j 

2 

8. 

36.9944 

6.0406 

13.8979 

''3j 

2 

9. 

41.0414 

6.7014 

0. 

''3j 

1 

10. 

44.7239 

7.3027 

36.5868 

''3j 

2 

11. 

47.0193 

7.6775 

10.8934 

3j 

2 

(a) 

We  used^® 

a(PbSe , 

T  -  300K)  - 

6.1243  A 

(b) 

Equation 

( D -  8 ) ,  and 

j  -  0,  ±  1, 

±2  ... 

NOTE:  The  allowed  g  values,  the  hexagon-tilt  angle  ♦pCg)  (see  Equation  3-1),  the  L  in 
units  of  2Ti/a*  and  the  hexagon  multiplieity  for  a  fee  lattiee  along  the  [111] 
direetion. 
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We  turn  now  to  the  determination  of  the  allowed  i  values  in  the  wave 
function  multipole  expansion,  Equation  (2-11)  ,  and  the  pseudopotential 
multipole  expansion  Equation  (2-9).  Since  the  Hamiltonian  is  invariant  under 
the  space  group  operators,  acting  on  ^piJ(?)  with  any  of  these  operators  yields 
another  energy-degenerate  wave  function  associated  with  another  Ic  vector.  The 
subgroup  of  operators  which  leave  a  Ic  vector  unchanged,  or  yield  an  equivalent 
vector  Tc  +  t  where  C  is  a  reciprocal  lattice  vector,  is  called  the  "group  of 
wave  vector  Ic" .  For  a  given  Ic,  the  energy-degenerate  wave  functions  'l'^j^(t) 
(a  =  1,2,...  is  a  degeneracy  index)  span  the  irreducible  representations  of 
the  group  of  the  wave  vector  Ic.  For  symmorphic  space  groups,  the  group  of 
wave  vector  Ic  coincide  with  a  subgroup  of  the  associated  point  group.  For 
non- symmorphic  groups  (i.e.,  groups  with  screw  axes  and/or  glide-reflection 
planes)  the  construction  of  the  group  of  wave  vector  Ic  is  more  complex,  in 
particular  for  Ic  vectors  on  the  Brillouin  zone  surface.  In  the  case  of  PbSe, 
the  space  group  0^  is  symmorphic  with  the  point  group  0^^. 

For  the  F-L  direction  in  the  Brillouin  zone  the  relevant  groups  of  wave 
vector  Ic  are  and  for  a  A  point  and  the  L  point,  respectively.  These 
groups  and  their  irreducible  representations  are  discussed  in  Appendix  A.  To 
demonstrate  the  determination  of  the  allowed  i.  values,  consider  the  two  one¬ 
dimensional  representations  at  a  A  point, (Table  3-2). 

Starting  from  a  general  multipole  expansion 


"Kr') 


'I'^(p,z')e 


(3-3) 


we  apply  and  Yq,  Equation  (A-1).  These  steps  give 
X  |'K?')>  =  X  I-  (p,z')  e^^^***  -  3"^ 


l'Kr')>  -  X  'I'p(P.z') 


(3-4) 


Yq  |'I'(?')>  -  I  ^^(p.z')  e’^^* 


e|«'(r')> 


I  ’l'n(p.z') 


(3-5a) 
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TABLE  3-2.  IRREDUCIBLE  REPRESENTATIONS  OF 


REPRESENTATION  -+  A-^  A2 

OPERATOR 


A3 


X 


0 


1 


1 


1  O' 

.0  1. 


1 


1 


■u)±^  0  ' 

.  0 


Y 


0 


1 


(Plane  Reflection) 


-1 


’0  1' 

.1  0. 


1  -1 


•  0 

0  , 


(a)  w  - 
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where 


e  -  +  1  for  , 

-  -  1  for  A2  .  (3-5b) 

•  »2ii 

+1,2 — 

Equating  both  sides  of  Equation  (3-4)  yields  e~  3  -  1  or 

j2-3m  ,  m-0,  +1,  +2,...  (3-6) 

The  same  procedure  when  applied  to  Equation  (3-5)  yields  the  i-parity: 

'I'_^(p,z')  -  e  *^(p.2')  (3-7) 


Combining  Equations  (3-6)  and  (3-7)  gives  the  multipole  expansion  of 
which  transforms  as  A^,  A2  under  the  operators; 


^'(r') 


1 1>  -  'I'q(p,z')  +  24'^(p,2')  cos  3(j)  +  24'g(p,z')  cos  64)  +  ... 


A^'  ^'(r')  -  1 1>  -  2i'I'2(p,z')  sin  34>  +  2i’l'g(p,z')  sin  6«^  +  ...  (3-8) 

The  above  exercise  demonstrates  how  the  symmetry  of  the  group  of  wave 
vector  Ic  severely  restricts  the  multipole  expansion  i  values  and  determines 
the  i-parity  (i.e.,  a  relation  between  and  ♦jj(pp,z')).  The 

treatment  of  the  two  dimensional  representation  A^  is  given  in  Appendix  B  with 
the  following  multipole  expansion 

|l>  -  ...  +  ,z' )e  +  ^^(p.z' )e^'^  +  4'^(p,z')e^^'^  +  ... 

A3: 

|2>  -  ...  +  'I'^(p,z  ’  )e'^^**’  +  'l'j^(p,z')e’^'^  +  P  •  ^  +  ...  (3-9) 


Note  that  both  in  Equations  (3-8)  and  (3-9)  the  p  dependence  is  retained  in  a 
multipole  'l'^(p,z').  When  the  p-dependent  Bessel  function  (see 
Equation  (2-lla))  is  factored  out,  the  i-parity  of  ♦^(g.z’)  are  those  listed 
in  Table  3-3. 

At  the  L  point  the  group  of  wave  vector  Ic  is  which  contains  in 

addition  to  the  operations  the  z' -inversion  (Equation  (A-1)). 

Consequently,  the  multipoles  have  in  addition  a  definite  z'-parity  (i.e.,  the 
phase  relation  between  'fjj(gp,z')  and  'l'jg(gp, -z ' ) )  .  Following  the  same  steps  as 
above  yields  the  sequences  of  allowed  i  values  for  all  representations,  see 
Table  3-4. 

The  pseudopotential.  Equation  (2-5),  is  by  construction  invariant  under 
all  space  group  operations,  in  particular  under  the  groups  of  wave  vector  Ic. 
Therefore,  it  transforms  as  the  scalar  representation  under  the  C^.^  and 
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TABLE  3-3.  MULTIPOLE  CONTENT  OF  THE  IRREDUCIBLE  REPRESENTATIONS 


i -parity 

IRREDUCIBLE 

REPRESENTATION  i-values^^^  d 


'^l 

ll>; 

i  -  3m 

(-1)^ 

1 

hi 

|1>: 

£=  3m,  m^O 

-(-1)^ 

1 

A3 

|1>: 

£  -  3m  +  1 

|2>: 

i  -  3m  -  1 

.(-i)“  oj 

2 

(a)  m  “  0 ,  +1,  +2,  .  .  .  . 


(b)  (’J'.^(g,z'))  s  n^(’If^(g,z'))  where  ('l'^(g,z')) 


♦^(g.z')  or 


’^^^hg.z'Y 


NOTE:  The  allowed  £  values  and  i -parity  for  the  group  C3V 
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TABLE  3-4.  MULTIPOLE  CONTENT  OF  THE  03^  IRREDUCIBLE  REPRESENTATIONS 


i-parity 


z- parity 


REPRESENTATION 

i 

-values^^^ 

n 

d 

Ll 

|1>: 

i  -  3m 

(-1)^ 

(-1)^ 

1 

L2 

|1>: 

i  -  3m,  ra^O 

-(-1)^ 

(-1)^ 

1 

/ 

Ll 

|1>: 

i  -  3m,  m?'0 

-(-1)^ 

-(-1)^ 

1 

l' 

|1>: 

i  -  3m 

(-1)^ 

1 

L3 

|1>: 

|2>: 

i  -  3m  +  1 

i  -  3m  -  1 

■0  (-1)“] 

.(-1)“  0, 

- 

■(-!)“  O' 

.0  (-1)“, 

2 

L3 

|1>: 

|2>: 

i  =  3m  +  1 

i  =  3m  -  1 

■0  (-!)“■ 

.(-1)“  0. 

■(-D®  O' 

P  (-1)®. 

2 

(a)  m  -  0,  +1,  +2 . 

(b)  (4'.^(g,z'))  =  n^d'^Cg.z'))  where  (’J'^{g,z'))  -  ^-^Cg.z')  or 

(c)  (^^(g.-z'))  =  )) . 


[«^2^(g,z') 


NOTE:  The  allowed  i  values,  i-parity  and  z-parity  for  the  group 
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operators,  i.e.,  as  and  respectively.  This  imposes  further 

restrictions  on  the  relevant  i -values,  as  demonstrated  now. 

Consider  for  example  the  pseudopotential.  For  our  working  example  it 
turns  out  that  at  most  the  four  lowest  g  values  contribute  (Section  5)  and 
only  0<i<6  need  be  considered  (Equation  (3-2a)).  On  the  other  hand,  the 
pseudopotential  transforms  as  the  representation,  hence,  according  to 
Table  3-4  i  -  0,  ±3,  ±6,  ...  .  Consequently,  the  independent  pseudopotential 
multipoles  have  only  two  £  values  e.g.,  £  -  0,3!  All  the  other  multipoles  are 
phase-related.  The  same  considerations  when  applied  to  the  wave-function 
(Section  5)  reduce  drastically  the  number  of  independent  multipoles. 
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SECTION  4 

THE  MULTIPOLES'  COUPLED  WAVE  EQUATIONS 


Having  established  the  symmetry  properties  of  the  multipole  expansion, 
our  goal  is  to  transcribe  the  Schrodinger  equation  for  into  equations 

for  the  multipoles  't^Cgp.z').  We  sketch  here  the  derivation  while  most 
details  are  given  in  Appendix  C. 

The  Schrodinger  equation  for  a  given  pseudopotential  v(r) ,  is 


r-2  ^'"O 
[v  +  ^ 


(En(k) 


v(r)) 


J  nk 


(r) 


(4-1) 


where  mQ  is  the  effective  electron  mass  (approximated  here  by  the  free 
electron  mass),  Eri(l^)  is  the  band  energy  for  a  given  band  index  n  and  momentum  ^ 
and  wave  function.  Since  is  a  scalar  under  all  three 

dimensional  rotations,  the  Schrodinger  Equation  (4-1)  in  the  primed 
coordinates  (Equation  (2-1))  takes  the  form: 

[  (v'^  +  (E^(k)  -  v'(?'))j  (?')  -  0  (4-2) 

Equation  (4-2)  is  now  transcribed  into  an  equation  for  the  cylindrical 
coordinates  multipoles  Equation  (2-2).  Expressing  the  la^lacian  in  the 
cylindrical  coordinates  and  introducing  the  multipole  expansion  of  the  product 


X(?')  -  v'(?')  'f;j!(r')  -  II  X^(gQ.z')J^(gQp)e 

’0 


ii4> 


we  obtain  after  some  projection  manipulations 


2  ^0  d^ 

4  dz  ‘ 


^"*0 

*^(go-z')  - 2  “  0 

ti 


(4-3) 


(4-4) 
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To  render  Equation  (4-4)  useful  it  is  necessary  to  express  X^(gQ,z')  in 
terms  of  the  underlying  multipoles  v^(gp,z')  and  'I'^(gp,z').  This  is  achieved 
by  means  of  a  decomposition  theorem  (Appendix  C)  which  states  that: 


e  T 


L^) 


3 

I 

m«-2 


■imS(G^) 


■i(i-m)6(l^) 


(4-5) 


In  Equation  (4-5)  and  are  two  transversal  reciprocal  vectors  of  lengths 

gp  and  gp,  respectively,  gQ.gp,  and  gp  add  up  to  the  triangle  defined  in 
Figure  4-1  and  A^  m^^O’^P’^F^  geometrical  coefficients  given  in  Table  4-1. 
In  terms  of  Equation  (4-5)  it  is  possible  to  express  X^(gQ,z')  Equation  (4-3) 
in  a  separable  form.  These  steps  yield  the  basic  wave  equations: 


)  + 


3  r  1 

^  ^  ^  ,  ^.m^SO-gp-gF)  [-  Vep-^'>J^.m(gF-^')  ’  » 

<7—  mss  -  ?  L  ^  J 


gp-Sp  m--2 


The  'I'^(gp,z')  satisfy  the  boundary  condition  implied  by  the  Bloch  theorem  in 
the  [111]  direction,  i.e.. 


«'.(g  ,z’  +  a  ) 


(4-7a) 


where  k',  and  a  are  given  by 


k' 


-kI-JT 

•a  ^ 


K 


3ti 


IT 


0  <  K  <  1 


(4-7b) 


where  the  L  and  a  A  point  correspond  to  K  -  1  and  K  0,  respectively. 


Equations  (4-6)  and  (4-7)  are  the  central  results  of  this  work.  These 
equations  are  exact  and  valid  in  the  [111]  direction.  The  symmetry  of  the 
bands  is  dictated  by  the  pertaining  i-values  sequence  (Tables  3-3,  3-4).  The 
participating  g-values,  in  ascending  order  (Table  3-1),  are  determined  by  the 
convergence  properties.  Equation  (4-6)  can  be  interpreted  as  a  set  of  coupled 
Kronig  Penney  models.^'  The  coupling  reflects  the  three-dimensional  nature  of 
the  Schrodinger  equation.  Note  that  the  diagonal  term  "-gQ^"  represents 
substraction  of  the  "transversal"  kinetic  energy  from  the  total  energy  E(lc)  of 
the  band.  This  is  expected  since  Equation  (4-6)  involves  only  z'.  The  manner 
in  which  Equation  (4-6)  is  solved  is  a  matter  of  convenience.  We  solved  it  by 
transcribing  it  into  a  secular  equation  in  a  plane  wave  basis  (Appendix  D) . 
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FIGURE  4-1.  THE  TRANSVERSAL  RECIPROCAL  VECTORS  PARALLELOGRAM 
DEFINING  THE  NOTATIONS  USED  IN  SECTION  4 
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TABLE  4-1.  THE  A-COEFFICIENTS 


gp.Bp 


TRIANGLE 

configuration 


Sp=gF=0 


A.AL.'-^'^  gp=gpi'0 


t  (-1)^ 


gp=^'-  &F-&0 


gp-BO'  6f“° 


A.AL /'''  A.AL,'^  '  gp^^O.  gp^^O  A„ 


(a)  c 


go  ■  sf  ■  Sp 

2  gp  gp 


I  c  I  =-  2  '  2  ’  ^  crystal  and  the  '111]  direction. 


.m^So’gp'SF^  “  i  B^^''’^exp(  imjarccos  (c)|  )j 

(((gp  +  gp  (t  -  i  ,  1-c^  ))/go)'^'  W(|c|)  -  j  for  |c|  =  l 


1  for  |c|<l 


(c)  A.AL,  =  Any  allowed  value. 


NOTE:  The  A - coef f ic lent s  The  triangle  configuration  corresponds  to  Figure  4-1  For  all 
case.s  not  quoted  in  the  Table  A^  ni(gO'6p’gF^  “ 
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SECTION  5 

BAND  STRUCTURE  OF  PbSe  IN  THE  F-L  DIRECTION 


In  this  section  we  test  the  multipole  expansion  convergence,  i.e.,  that 
a  small  number  of  multipoles  yields  convergent  results  for  Ej^(lc) .  The 
quotations  added  to  "PbSe"  are  a  reminder  that  the  important  spin-orbit 
interaction  is  not  included  at  this  stage.  To  keep  the  discussion  focused, 
the  section  is  subdivided  into  subsections,  each  devoted  to  a  specific  aspect 
of  the  calculation.  Further  details  are  given  in  Appendix  D. 


THE  PSEUDOPOTENTIAL'S  MULTIPOLES 


The  input  to  the  wave  equations  (Equation  (4-6))  are  the 
pseudopotential's  multipoles.  The  empirical  pseudopotential  for  PbSe  is 
comprised  of  v(2)  such  that  |2|^<  12  (2ii/a)^  (Equation  (2-4))  implies  that 
59  such  vectors .contribute .  When  transformed  to  the  primed  coordinates  the 

the  first  four  g-hexagons,  i.e.,  gp  -  0,  1.67,  2.90, 

The  discussion  in  Section  3  implies  that  only  i-0,3  need 

be  considered. 


^  MW  ill 

corresponding  span 
3.35  A'^,  Table  3-1. 


The  Vg  and  v^  multipoles  are  depicted  in  Figures  5-1  and  5-2, 
respectively.  Note  that  V3(g(3),z')  -  0.  This  is  a  consequence  of  the 
following  juxtaposition:  Table  3-4  implies  in  this  case  that 
v_2(g(3),z')  "=  — V2(g(3) ,z ' ) .  On  the  other  hand,  Equation  (3-2a)  with 
^0^S(3))  =  n/6  (Table  1)  gives  v,3(g(3),z')  -  V3(g(3),z').  Hence 
V3(g(3),z')  =  0.  Figure  5-3  displays  the  full  pseudopotential  along  the  [111] 
direction,  i.e.,  (Equation  (2-9a)): 

v'(x'  =  y'  -  0,  z')  -  X  VQ(gp.z')  (5-1) 

gp 

The  deeper  well  at  z '  -  a  /2  indicates  that  it  is  the  anion  location  since 
electrons  in  the  lowest  bands  accumulate  there.  Correspondingly,  the  Pb  atom 
is  at  the  origin.  This  outcome  is  also  inferred  from  the  values  of  v(3) ;  see 
the  remark  subsequent  to  Equation  (D-6). 

The  magnitude  of  the  multipoles  v^(gp,z')  is  indicative  of  the  important 
gp  values.  This  is  evident  from  Equation  (4-6),  where  these  multipoles 
provide  the  coupling  between  the  'I'^(gp,z').  As  Figures  5-1  and  5-2  show,  the 
gp  =  g(2)  multipole  is  particularly  large.  Hence,  this  g-value  is  definitely 
needed.  The  g  -  g(l)  must  be  included  since  it  enters  when  two  equal  g-values 
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FIGURE  5-1.  THE  I  =  0  MULTIPOLES  OF  THE  PbSe  PSEUDOPOTENIAL 
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FIGURE  5-2.  THE  #  =  3  MULTIPOLES  OF  THE  PbSe  PSEUDOPOTENIAL 
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FIGURE  5-3.  THE  PSEUDOPOTENTIAL  ALONG  THE  z’flU]  AXIS  (EQUATION  (5-1)) 
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raultipoles  combine.  The  g  -  g(3),  g{4)  are  border  cases  (see  Figure  5-1)  and 
their  significance  can  be  assessed  numerically.  All  in  all,  it  is  necessary 
to  check  four  g-values,  i.e.,  g(l),  g(2),  g(3),  g(4). 


CONVERGENCE 

The  usefulness  of  the  wave  equations  (Equation  (4-6))  hinges  on  the 
convergence  properties  of  the  multipole  expansion.  We  examine  this  issue  by 
calculating  the  bands  along  the  entire  F-L  direction  in  the  Brillouin  zone. 

The  relevant  parameter  domain  is  as  follows.  As  argued  in  subsection 
Equation  (5-1),  the  potentially  important  g  values  are  the  lowest  four 
(Table  3-1).  As  Tables  3-3  and  show  (for  these  g-values),  only  i  -  0,3 
need  be  considered  fqr  the  ,  L2  and  Aj^,  representations .  By  the  same 

token,  for  the  L2 ,  and  A2  representations,  only  i  -  -2,  +1  need  be 
considered . 

In  evaluating  the  calculated  bands,  we  look  for  a  good  match  with  the 
known  band  gap  at  the  L-point,  the  correct  number  and  sequence  of  the  bands 
(as  obtained  in  other  band  calculations),  and  the  correct  band  degeneracies  at 
the  F-point.  The  calculations  were  (Jone  considering  g-values  in  ascending 
order.  The  calculation  of  the  ,  L2  bands  and  L^.L^'  bands  decouple  and, 
hence,  they  are  considered  separately.  The  results  are  summarized  in 
Table  5-1  for  the  L-point  and  in  Figure  5-4  for  the  entire  F-L  line. 

For  the  L^ ,  L2  bands  the  lowest  two  g-values  at  least  are  needed,  i.e., 
g  -  g(l),g(2)  (Table  3-1).  This  three-multipole  base  yields  a  band  gap 
AE-.921  eV,  too  large  in  comparison  with  the  experimental  value 
AE(PbSe,  T  -  4.2K)  -  .165  eV.^^  as  Table  5-1  shows,  adding  the  (i-0,  g-g(3)) 
multipole  reduces  the  gap  to  AE-.531  eV,  i.e.,  this  multipole  is  important. 

On  the  other  hand,  the  (i-3,  g-g(3))  multipole  does  not  change  the 
eigenenergies ,  i.e.,  it  is  decoupled  from  all  other  multipoles  and,  hence, 
ignored  hereafter.  Adding  the  g“g(4)  multipole  reduces  AE  further  by  .1  eV. 
This  indicates  convergence.  The  convergent  result  for  the  gap,  however,  is 
larger  than  the  experimental  value.  This  outcome  is  expected  since  the  spin 
orbit  interaction,  known  to  reduce  the  gap,  is  omitted  here.^^  For  the 
(L^,  L^)  representations,  a  similar  pattern  emerges  (Table  5-1).  Including 
only  the  g=g(2)  multipoles  is  insufficient.  The  g-g(3)  multipoles  give 
substantial  corrections  to  the  pertinent  gap,  and  the  g-g(4)  multipoles  give  a 
relatively  smaller  correction,  indicating  convergence.  We  did  not  check  the 
effects  of  modifying  mg,  the  effective  mass.  The  number  and  sequence  of  the 
calculated  bands  is  the  correct  one.^^ 

The  expansion  of  the  above  analysis  to  the  entire  F-L  line  yields 
Figure  5-4.  Three  representative  cases  are  displayed.  Figure  5 -4a ^ shows  the 
results  of  a  minimum-base.  This  includes  g-g(p,g(2)  for  the  Lj^ ,  L2 
representations  and  g-g(2),g(3)  for  the  (L^,  Lj)  representations.  The  results 
reproduces  the  correct  trends  and  the  degeneracies  at  the  F-point.  However, 
the  calculated  band  gap  (AE-.921eV,  Table  5-1)  is  much  larger  than  the 
experimental  value.  This  implies  that  more  multipoles  are  required  to 
adequately  describe  the  L^ ,  L2  representations. 
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TABLE  5-1.  CONVERGENCE  OF  BAND  ENERGIES  IN  PbSe  AT  THE  L-POINT(a) 
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NOTE:  The  six  lowest  bands  of  PbSe  at  the  L-point.  The  upper  part  of  the  table  lists  the  total  number  of 
imiltipoles  n  and  the  (#,g)  pairs  considered.  The  g-values  are  denoted  as  the  appear  in  Table  3-1,  e.g., 
g(l)  is  denoted  as  1  in  the  g-column.  An  asterisk  marks  a  multipole  included  in  the  calculation.  The 
lower  part  of  the  table  lists  the  corresponding  symmetries,  band  energies,  and  the  bad  gap  AEg  at  the 
Fermi  level. 
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FIGURE  5-4.  THE  CALCULATED  PbSe  BAND  STRUCTURE  IN  THE  T-L  DIRECTION 
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In  Figure  5-4b  we  display  a  "compromise"  base;  the  trends  and  bands 
arrangements  are  as  in  Figure  5-4.  In  addition,  the  calculated  L-point  band 
gap  is  substantially  reduced  in  comparison  to  Figure  however,  the  F- 

point  degeneracies  are  not  reproduced  since  the  bands  are  not  evaluated 

to  a  similar  degree  of  accuracy.  Figure  5-4c  shows  the  convergent  calculation 
with  a  calculated  gap  of  AE-.293eV,  (Table  5)  and  the  corect  F-point 
degeneracies.  Six  multipoles  are  included  in  this  calculation  for  both  L2 
and  ,  L-j  representations.  This  result  demonstrate  that  good  convergence  is 
achieved  for  the  entire  F-L  line  with  a  small  number  of  multipoles. 
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SECTION  6 

SUMMARY  AND  DISCUSSION 


We  developed  a  method  to  calculate  bands  along  a  high-sjnnmetry  axis 
based  on  the  usage  of  cylindrical  coordinates.  This  representation  leads  to 
the  general  expansion  of  a  lattice -periodic  function  in  terms  of  single - 
variable  cylindrical  multipoles. 

The  two  central  results  of  this  work  are  (1)  the  derivation  of  the  wave 
equations  set  for  the  multipoles  and  (2)  the  demonstration  of  good  convergence 
of  the  multipole  expansion  in  a  test  case.  Using  sjrmmetry  considerations  the 
three-dimensional  Sc'^  r  inger  equation  with  a  locc  1  potential  is  exactly 
reduced  to  a  set  of  coupled  ordinary  differential  equations  for  the 
cylindrical  mult...p  ...es.  This  reduction  employs  new  A-coef f icients  which 
express  the  constraints  of  the  transversal  (two-dimensional)  vector  addition, 
in  analogy  tn  the  Clebsh-Gordan  coefficients  in  the  context  of  three 
dimensional  angular  momentum. These  coefficients  vary  with  the  point  group 
and  high  symmetry  axis.  The  expressions  for  the  [111]  direction  in  a  fee 
lattice  are  derived. 

We  have  also  demonstrated  good  convergence  of  the  cylindrical  multipole 
expansion  in  a  simple  case,  i.e.,  PbSe  along  the  [111]  direction.  These 
multipoles  depend  on  two  discrete  indices,  denoted  by  i  and  g.  The  former 
plays  the  role  of  projected  (longitudinal)  angular  momentum.  The  latter  is  a 
transversal  momentum.  Group  theoretic  considerations  limit  the  number  of 
distinct  i -values  to  two .  thus  rendering  the  multipole  expansion  in  effect  to 
an  expansion  in  g.  Gonsequently ,  the  convergence  is  determined  by  how  many  g 
values  contribute.  As  argued  in  the  INTRODUCTION,  the  convergence 
considerations  ^  not  depend  on  the  z-extension  of  the  unit  cell.  Hence,  the 
present  approach  is  expected  to  be  particularly  advantageous  for  highly 
anisotropic  structures,  e.g.,  superlattices.  This  conclusion  can  also  be 
stated  as  follows:  Whereas  the  number  of  basis  functions  in  a  bulk-method 
grow  as  N  ,  where  N  is  the  number  of  basis  functions  (e.g.,  plane  waves) 
needed  for  each  direction  in  space,  the  number  of  basis  functions  in  present 
method  grows  as  mN,  where  m  is  a  small  number  (<6). 

The  present  method  bears  some  similarity  to  the  method  of  symmetrised 
plane  wave  method  (SPWM).^^  In  both  methods  the  symmetry  at  the  fe-point  is 
utilized  to  construct  "multipoles"  to  reduce  the  dimensionality  of  the  secular 
equation.  The  present  approach  has  the  following  new  features.  It  goes  one 
step  further  by  introducing  the  cylindrical  coordinates  system.  As  a  result, 
the  present  multipoles  depend  on  one  coordinate  while  the  SPWM  multipoles 
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depend  on  three  coordinates.  Secondly,  we  derived  the  set  of  coupled  wave 
equations,  Equation  (4-6),  for  the  multipoles.  These  are  exact  and  do  not 
depend  on  a  specific  basis  (e.g.,  plane  waves).  Future  work  may  employ  more 
efficient  algorithms  than  the  one  employed  in  Appendix  D.  Third,  the 
convergence  properties  do  not  depend  on  the  extension  of  the  unit  cell  in  the 
high-symmetry  direction.  On  the  other  hand,  both  SPWM  and  the  present  method 
cannot  be  readily  used  for  first  principle  calculations  since  they  do  not 
yield  the  bands  throughout  the  Brillouin  zone.  Thus,  application  of  the 
present  method  is  primarily  expected  in  phenomenological  calculations. 

The  present  approach  suggests  extensions  in  several  directions.  The 
inclusion  of  the  spin  orbit  interaction  is  one.  Other  possibilities  are 
extensions  to  different  high-sjrmmetry  axes,  the  band  structure  at  the  F  point, 
and  different  point  groups.  We  hope  to  address  these  matters  in  future 
publications . 
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APPENDIX  A 

IRREDUCIBLE  REPRESENTATIONS  OF  THE  D3J  and  03^  GROUPS 


For  our  purpose  it  is  particularly  convenient  to  use  the  cylindrical 
realizations  of  €3^  and  03^  groups 


£3v 

Xq  'Kcj))  =  'Ktt)) 

X+i  ±  2ii/3) 

Yq 

'K4))  =  ^K-cl)  +  2ii/3) 


^3d 

Xq 

±  it/3.  -z) 
X+2  ’P(4>.z)  =  ±  2ti/3,  z) 

X3  'i'((|>,z)  =  +  n,  -z) 

Xq  'I'(4>,z)  -  ’K-cj).  z) 

\l  «'(<!>, z)  -  ±  it/3,  -z) 

Y^2  'I'(4'.z)  “  ±  2ti/3,  z) 

X3  ^'(4),z)  -  ^(-<|)  +  n,  -z) 


(A-1) 


In  Equation  (A-1)  ^^(<1))  and  5'(4>,  z)  are  arbitrary  functions  and  (}>  is  the 
cylindrical  azimuthal  angle.  Note  that  Yq  is  a  reflection  about  a  plane  which 
includes  the  z-axis  and  X3  is  the  inversion  about  the  origin. 

To  derive  expeditiously  the  irreducible  representations  (F),  it  is 
necessary  to  find  a  convenient  basis  set  Ujjj.  Once  such  a  set  is  selected,  (F) 
is  determined  according  to^'^ 

A  A 

R-  u  -  Z  Uu(r(R- ))v  (A-2) 

1  m  ,  k  1  km 
k 
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A 

where  is  a  group  operator,  Equation  (A-1).  For  the  group  there  are 
three  irreducible  representations  for  which  we  choose  the  following  bases: 


All  ll>  -  e^i*  +  e-3i«l> 

A2:  |1>  -  e3i4>  . 

A3:  |l>  =  ei'»» 

|2>  - 


for  €3^ 


The  ensuing  irreducible  representations  are  given  in  Table  A-l. 


(A-3) 


For  the  group  there  are  six  irreducible  representations  for  which 

the  following  bases  have  been  chosen: 


l{:  |1>  =  (e^^***  -  U(z), 

L2;  |1>  -  (e^^***  +  U(z), 

Lj^:  |1>  -=  U(z), 

L2:  |1>  -  U(z). 

L3:  |1>  -  e^***  U(z) 

|2>  -  e"^'*’  U(z)  ,  U(-z 

L3:  |1>  =  e^***  U(z) 

|2>  -  e'^***  U(z)  ,  U(- 


U(-z)  =  U(z) 

U(-z)  -  U(z) 

U(-z)  -  -U(z) 

U(-z)  -  -U(z) 

-  U(z) 

,  for  D3jJ 

)  -  -U(z)  (A-4) 


The  corresponding  irreducible  representations  are  given  in  Table  A-2. 
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TABLE  A-1.  IRREDUCIBLE  REPRESENTATIONS  OF 


REPRESENTATION  -»  A;i^ 

OPERATOR 


X 


0 


1 


A3 


'1 

0 


1 


1 


Yq  1  -1 

(Plane  Reflection) 


’0 

1 


1  -1 


(a)  u,  -  e2^i/3 
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TABLE  A-2.  irreducible  REPRESENTATIONS  OF  D 


REPRESENTATION  -►  L^^  L2  L2  Lj^  L3 

OPERATOR 


(Plane  reflection) 
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(p,z')  -  (B-5) 

For  the  group  we  use  Table  A- 2  and  apply  the  operators  X^l.  X3  and 
Yq  (Equation  (A-1)).  This  gives  for  the  one  dimensional  representations 

X+i  ll>  =  I  ’I',  e  |1>  (B-6) 

Jl^-ao 

where  (Table  A- 2) 


e  =  1  for  L^, 

=  -1  for  L^, 


(B-7) 


In  a  similar  way 


X3  |1>  “  I  'I'  /p.-z')e 

i— « 


ii  (<|)+n) 


1> 


(B-8) 


and 

00 

Yq  |1>  '  Z  «  (p,z')e'^^‘*’  -  6  |1>  (B-9) 

where 

6  ■=  1  for  L^, 

=  -1  for  L^,  L2  (B-10) 

Equations  (B-6)  and  (B-7)  give 

'I'^(p.-z')  -  e  ♦^(p.z')  (B-11) 

The  LHS  of  Equation  (B-11)  holds  for  i  -  3m,  where  m  -  0,  +1,  +2 . 

The  RHS  of  Equation  (B-11)  gives  the  z-parity  rule,  i.e.,  a  relation  between 
'J'^(p,-z')  and  'I'^(p,z').  Finally,  Equation  (B-9)  yields  the  f-parity  rule 

^(p,z')  -  6  ♦^(p,z')  (B-12) 

Note  that  Tables  3-3  and  3-4  of  the  text  refer  to  4’^(g,z'). 
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The  two-dimensional  representations  Lo,  are  treated  as  above. 
Applying  ,  X-j  and  Yq  yields  equations  sufficient  to  determine  the  Z  values 
and  parities. 
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APPENDIX  C 

DERIVATION  OF  EQUATION  (4-6)  AND  THE  DECOMPOSITION  THEOREM  EQUATION  (4-5) 


To  arrive  from  Equation  (4-4)  of  the  main  text  to  the  multipoles  wave 
equation  Equation  (4-6),  it  is  necessary  to  express  X^(gQ,2')  in  a  separable 
form.  The  explicit  expression  of  Xjg(g^,z' )  (Equation  (4-3))  is; 

X^(gQ,z')  =  i^  X  v'(G')  <K'(L')  ^T^  (C-1) 


|Gt+  l,j,|  -  gQ 


where  'I'V(L')  =  'I''(L')  and 
nk 


G'+  L' 

X  X 


cos  6(G,^+  Lp 


,  sin  &(G^+  Lp 


G'+  L' 

y  y 


(G'+  L')^+  (G'  +  L')^ 

XX  y  y 


(0-2) 


(for  gQ  =  0  only  i  -  0  is  allowed  and  6(2^'+  tj')  -  0).  The  constraint  on  the 
summations  in  Equation  (C-1)  reflects  the  geometry  of  a  triangle  formed  by 
vector  addition  of  2^'  and  TLj' ,  see  Figure  4-1  in  the  main  text. 

The  constraint  on  the  summations  in  Equation  (C-1)  can  be  written  in  a 
more  convenient  form  as; 


G'  ,L' 


-  I  I 

Sp’Sp  G' 


(C-3) 


Gt+  -  gQ 


|G^I  -  gp  |L^I  -  gf 

|g;  +  L'l  -  gQ 


In  Equation  ((|:-3)^the  sum  over  C',  2'  is  converted  ^o  a  suip  over  all  possible 
hexagons  of  2^,  Cj.  However,  not  all  hexagons  of  llj.  and  2.j.  are  mutually 
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compatible:  Since  Up,  and  Ilp  +  2^.  form  a  triangle,  only  those  gp  and  gp 

which  form  a  triangle  with  the  third  side  of  length  gQ  are  ^ccep^able.  The 
other  gp,  gp  combinations  are  excluded  by  the  constraint  [Up  +  2p|  “  gg-  This 
constraint  is  quantified  by  considering  the  triangle  in  Figure  4-1.  First  we 
consider  non-degenerate  triangles  (i.e.,  when  gg ,  gp,  gp  0  and  no  angle  in 
the  triangle  is  0  or  n).  All  other  cases  are  treated  thereafter. 

The  cosine  theorem  for  the  triangle  in  Figure  4-1  gives 

2  2  2 
gg  ■  Sp  *  gp 

cos  [6(Gp  -  6(L,^)] - - c  (gQ.gp.gp)  =  c  (C-4) 

Hence  a  triad  gg.gp.gp  represents  a  triangle  if,  and  only  if,  |c|<l.  The  RHS 
of  Equation  (C-3)  then  takes  the  form 

I  I  I  (C-5) 

gp.gp  g<  L' 

1  (  gg  >  gp  >  gp  1 '^1  gp  l^-pl  “  gp 

where  in  Equation  (C-5)  indicates  a  summation  only  on  gp,gp  such  that  for  a 
given  gg  the  inequality  |c|<l  holds. 

The  other  source  of  "coupling"  in  Equation  (C-1)  is  in  the  shift  phase 
6(2,p  +  .  To  simplify  this  term  we  derive  first  an  identity.  Projecting 

the  triangle  in  Figure  4-1  on  the  x'  and  y'  axes  and  adding  the  two  equalities 
with  the  imaginary  unit  factor  gives  the  identity 


^±i6(G,|,+  L^) 

'  So 


^±i6(G,^) 


(C-6) 


On  the  other  hand,  inversion  of  Equation  (C-4)  gives 

6(G.^)  -  6(L.j,)  -  ±  a  ,  a  -  jarccos  (c)|  (C-7) 

where  the  two  signs  in  Equation  (^-12)  correspond  tq  the, two  possible 
triangles  (interchanging  2^  and  C^)  with  the  same  2p  +  tj.  Combining 
Equation  (C-4)  and  (C-6)  with  the  minus  sign  gives  the  identity: 

-i[6(G;  +  l;)  -  6(l;)] 

e  1  T  i 

-  +  b('^6  (C-8) 

S(Gj)  -  6(L,^)  -  a,0  6(G,^)  -  6(L,^)  +  a,0 
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where 


^  [gp  +  gp(c  -  i  Jl-c2  )],  (C-9) 

The  Kronecker  6  are  inserted  in  Equation  (C-8)  since ^the  LH^  is  under  the 
sununation  in  Equation  (C-1)  and  the  shift  angles  6(2^),  6(nj)  are  discrete. 
Furthermore,  since  on^O  for  a  non-degenerate  triai;igle  ^t  follqws  that  only  one 
term  in  Equation  (C-8)  contributes  for  given  6(2^.  +  tl^.) ,  .  Consequently 

raising  both  sides  of  Equation  (C-8)  to  the  i-th  power  one  obtains  the  key 
result : 


^-if  [6(G,J,  +  L,^)  -  6(L^)1 


6 

6(G,^)  -  6(L,^)  -  a.O 


& 

&(G,j.)  -  &(t|,)  +  a,0 


(C-10) 


where 


(b5^>)* 


(C-11) 


To  express  the  Kronecker  6  factors  in  Equation  (C-lO)^in  a  separable 
form,  consider  their  arguments.  Since  the  shift  angles  6(2^),  &(tj)  obey 
Equation  (3-1)  we  have 


6(G,^)  -  6(L,|,)  ±  a  -  n/3  (j-n)  +  4>o(gp)  *  “ 

=  n/3  (j-n)  +  A-  (C-12) 

■A* 

Now  Equation  (C-10)  is  under  the  Gj,  Lj.  suipmationf  of  Equation  (C-5).  This 
guarantee  that  at  least  one  triangle  Cj,  and  +  tj,  is  formed  with 

+  C.j.|  -=  gQ.  Therefore,  as  function  of  j  and  n,  the  RHS  of  Equation  (C-12) 
will  vanish  for  certain  j  and  n.  Hence  A~  must  equal  a  multiple  of  n/3  and 
the  LHS  of  Equation  Equation  (D-12)  is  always  a  multiple  of  n/3.  Obviously 
this  multiple  is  <6.  Therefore  the  Kronecker  6  in  (C-10)  can  be  expressed 
using  the  identity^^ 


,3  .n 

6  -  7  I  e'h  .  |s|<6  (C-13) 

^  s,0  m=-2 

Replacing  n/3  s  in  Equation  (C-13)  by  the  LHS  of  Equation  (C-12),  the 
Kronecker  &  factors  in  Equation  (C-10)  attain  a  separable  form  and  the 
decomposition  theorem  Equation  (4-5)  follows.  Expressing  Equation  (C-10)  in 
this  form  and  inserting  into  Equation  (C-1)  and  (C-5)  yields  the  wave  equation 
Equation  (4-6)  . 


C-3 


NAVSWC  TR  91-324 


A  similar  treatment  applies  for  the  degenerate- triangle  cases 
(Figure  4-1).  When  gQ  -  0,  only  i  -  0  is  allowed  and  6(G.^  +  L.^)  -  0  in 

Equation  (C-1).  Furthermore,  since  |G,|,  +  L.|,|  -  0,  it  follows  that 

|G.^|  -  |L.^|  -  g  and  the  vectors  G.^,,  L.^  point  back  to  back.  Therefore, 


Equation  (C-1)  takes  the  simple  form: 


X  (g  -  0,  z')  -  I  v'(G')  4<'(L')  e^^^z  ^z 

G'  ,  L' 


z  z 


+  I 

gp'O 


I  I 

G'  L' 

|2;|  =  g  |i:;i  =  g 


v'(G')  ^r'(L')  e^^^z  ^z  * 


+  6 


(C-14) 


6(G.j,) -6(1^) -71,0  6(G,|.)-6(L.^)+71,0 


The  first  term  in  Equation  (C-14)  corresponds  to  a  triangle  degenerated  into  a 

point,  i.e.,  -  0  -  g.  The  subsequent  lines  corresponds  to  a  triangle 

collapsed  into  a  segment.  As  can  be  checked  by  counting,  it  is  necessary  to 
include  the  two  Kronecker  6  terras.  However,  both  Kronecker  6  factors  are 
represented  in  this  case  by  the  same  expansion  Equation  (C-13).  This  is  the 
origin  of  the  "1/6"  in  Table  4-1  of  the  main  text. 

The  cases  gp“0,  gp“go  gp“g0’  Sp”®  trivial:  Obviously  in  these 

cases  S(G.j,  +  L.j,)  -  &(L.^)  and  6(G.^) ,  respectively  and  Equation  (C-1)  is  already 

in  a  separable  form.  For  the  other  case  with  c  -  -1,  two  Kronecker  6  factors 
contribute  in  Equation  (C-10);  however,  again  both  are  represented  by  the  same 

separable  expansion  Equation  (C-13).  For  the  case  with  c-1,  G,^  and  are 

collinear  and  point  at  the  same  direction.  Hence,  Equation  (C-10)  involves 
only  one  Kronecker  6,  which  leads  to  the  W(|c|)  -  1/2  factor  in  Table  4-1 
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APPENDIX  D 

NUMERICAL  SOLUTION  OF  THE  WAVE  EQUATIONS  (4-6) 


THE  PSEUDOPOTENTIAL  v(r) 

For  a  crystal  with  atoms  per  unit  cell,  each  of  specie  a,  the 
lattice-periodic  pseudopotential  is  given  by 


v(r)  =  X  (D-1) 

^  ,  a  a 

a  -=  1 

where  is  the  location  of  atom  a  with  respect  to  a  chosen  origin  in  the  unit 
cell  and  v^(t)  is  the  "elemental"  pseudopotential.  Since  Equation  (D-1)  is 
lattice-periodic,  it  has  the  (truncated)  fourier  expansion 


n. 


v(r)  =  I'  e 
-♦ 

G 


,  iG*r 


I  S^<G)v„(0) 

0(»1 


(D-2) 


where  the  structure  factors  Sq^(2)  and  the  "elemental"  parameters  v^(2)  are 
defined  by:^ 


n  -f 

S  (G) - e  a,  V  (G)  -  —  fdr  e  v  (r)  (D-3) 

an.  a  n  ^  a 

A  n 

In  Equation  (D-2)  2  has  a  length  smaller  that  a  given  cutoff,  and  in 
Equation  (D-3)  fJ  is  the  compound's  unit  cell  volume.  The  origin  of  the  n^ 
normalization  factor  in  Equation  (D-3)  is  evident  by  considering  the  n^  -  2 
case:  For  the  elemental  v^(^) .  the  corresponding  volume  factor  in 

Equation  (D-3)  is  that  of  the  elemental  solid.  However,  for  two  atoms  per 
unit  cell,  n  is  about  twice  the  volume  of  the  elemental  unit  cell.  Therefore, 
to  associate  v^(2)  with  the  elemental  parameters  it  is  necessary  to  introduce 
the  n^  factors  as  in  Equation  (D-3). 

For  the  special  case  of  n^^  -  2,  Equation  (D-2)  is  rewritten  in  the  form 
2 

I  So^(G)v^(G)  -  Sg(G)Vg(G)  +  S^(G)v^(G)  (D-4) 

a  =  1 
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where 

VgCG)  =  I  [Vj^(G)  +  V2(G)]  ,  v^(G)  -  i  [vj^(G)  -  v^CG)] 

Sg(G)  =  S^(G)  +  S2(G)  .  S^(G)  =  S^(G)  -  S^(G)  (D-5) 

For  the  rocksalt  structure  with  two  atoms  per  unit  cell  =  (0,0,0), 

T2  =  a( 1/2 , 1/2 , 1/2 ) ,  and  2  is  given  in  Equation  (2-A).  Consequently  the 
structure  functions  Equation  (D-5)  take  the  very  simple  form: 

85(0)  =  I  (1  +  exp[-in(n^+  +  n^)] 

S^(G)  -  I  (1  -  exp(-in(nj^+  n2  +  n^)]  (D-6) 

The  values  v^(2) ,  Vg(2)  for  the  corresponding  non  vanishing  structure 
functions  Equation  (D-6)  are  fitted. These  are  the  values  used  hereafter. 
By  comparing  the  quoted  values  of  Vg,v^  with  the  "elemental"  values  v^  for  Pb 
and  Se^*^  we  conclude  that  Pb  is  at  the  origin  and  Se  is  at  the  midpoint  of 
the  [111]  diagonal  (see  also  Section  3  and  Figure  5-1  of  the  main  text). 


NUMERICAL  DETAILS 

A  convenient  method  for  solving  the  wave  equations,  Equation  (4-6),  for 
the  eigenvalues  and  eigenfunctions  is  transforming  them  to  a  secular  equation 
in  a  plane  wave  representation.  Direct  numerical  integration  fails  for  the 
very  low  bands  where  the  equations  become  too  stiff  for  the  standard 
integration  routines. 

From  the  form  of  the  multipole  expansion  Equation  (2-12)  (main  text)  it 
follows  that 


ik'z' 


b(g,L^)  e 


iL'z' 

z 


(D-7) 


Martinez,  G.,  Schluter,  M. ,  and  Cohen,  M.  L. ,  Phvs .  Rev.  B.  Vol .  11,  No. 
#,  1975,  p.  651. 

Kohn,  S.  E. ,  Yu,  P.  Y. ,  Fetroff,  Y.  R. ,  Shen,  Y.  Tsang,  and  Cohen, 

M.  L.  ,  Phvs .  Rev .  B .  Vol.  8,  No.  y/,  1973,  p.  1477. 
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$ 

where  b(g,L2)  are  coefficients  to  be  determined.  Obvious!',  ,  Equ^.^  ion  (D-7) 
has  the  correct  Bloch  periodicity.  Furthermore,  invoking  Equations  (2-3),  and 
(2-4)  yields 


L '  =  r*  N  ,  N  =  n-  +  n„  +  r,  (D-8) 

z  a’^  i  3 

where  a*  =  a  J  3  ,  Equation  (4-7b)  ,  and  n^^,  n2,  n^  =  0,  +1,  +2,  ...  from 

Equation  (2-4).  Equation  (D-8)  implies  that  for  a  given  g  the  sequence  of 
allowed  N(n)  values  is  restricted  to  the  series  given  in  Table  1-1  of  main 
text.  These  values  are  either  of  the  form  3m  (m  is  an  integer)  or  the 
complementary  set  of  integers.  Imposing  these  restrictions  on  N(n)  is 
essential  to  avoid  spurious  degeneracy  between  the  bands  for  K  =  1  and  K  =  1/3 
and  to  reduce  the  size  of  the  secular  equation.  Hence,  the  plane  wave 
expansion  used  is: 


N 

w 

=  I  ^n^®^  (N(n)  -I-  I  K)]  (D-9) 

n=-N 

w 


and  N(n)  is  specified  for  each  g-value  in  Table  1-1. 

Inserting  Equation  (D-9)  into  the  wave  equations  and  projecting  the 
individual  plane  waves  yields  a  secular  equation  for  matrix 

elements  <n|V|^(gp,z ' ) |n'>  are  evaluated  using  a  Fast  Fourier  Transform 
algorithm.  Tha  eigenvalues  are  obtained  by  searching  for  E^(^)  which  make  the 
secular  equation  determinant  vanish.  Unphysical  Ej^(k)  can  be  identified  by 
considering  the  coefficients  b^(g).  For  an  unphysical  E^(lc)  these 
coefficients  nearly  vanish  for  one  or  more  of  the  multipoles  or  they  do  not 
vanish  for  the  proper  N(n) . 


At  the  L-point  (K  =  1,  Equation  (3-7b),  main  text)  every  multipole  has  a 
definite  z' -parity.  Table  3-4.  Combining  this  with  the  Bloch  periodicity 
Equation  (4-7)  (main  text)  yields 

'if^(gp,  z'  =  a*/2-C)  =  -^2  ^i^SF’  “  a*/2-t-C)  (D-lOa) 

where 


-a*/2  <  C  <  aV2  (D-lOb) 

Equation  (D-10)  implies  sign-relations  between  the  proper  pairs  of  b^(g) •  Iri 
this  manner  we  distinguish  between  the  Lp  and  bands.  For  K  ^  1,  all  bands 
belong  to  the  Ap  symmetry  (compare  Tables  3-3  and  3-4  (main  text)  to  obtain 
the  compatability  relations  from  the  i-parity).  Note  from  Tables  3-3,  3-4 
that  including  the  f=0  multipole  is  essential  to  differentiate  between  the  Lp, 
L2  and  the  L| ,  L2  representations. 

Our  experience  indicates  that  taking  -  7  in  Equation  (D-9)  is 
sufficient . 
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